/********************************************************************
 *
 *  psample_base.h
 *
 *  Some basic devices to support the implementation of probabilistic
 *  sampling.
 *
 *  Created by Dahua Lin, on Nov 6, 2010
 *
 ********************************************************************/

#ifndef DDP_PSAMPLE_BASE_H
#define DDP_PSAMPLE_BASE_H

#include <cmath>

namespace ddp
{


/**
 * get the index of interval that a value is in
 *
 * @param n         the number of consecutive intervals
 * @param redges    the right edges of intervals
 * @param v         the test value
 *
 * @return the index of the interval that v is in
 */    
inline int get_interval(int n, const double *redges, double v)
{
    int i = 0;
    while (i < n && v > redges[i]) ++i;
    return i;
}


/**
 * compute cumulative sum
 *
 * @param n     the number of elements in a
 * @param a     the input array 
 * @param cs    the output array
 *
 * @return the sum of elements in a
 */
inline double cumsum(int n, const double *a, double *cs)
{
    double s = 0;
    for (int i = 0; i < n; ++i) cs[i] = (s += a[i]);
    return s;
}


/**
 * draw a sample from a discrete distribution
 *
 * @param n     the number of classes
 * @param p     the probabilities
 * @param u     a uniformly distributed value
 * @param F     the cumulative distribution (output)
 *
 * @return the sampled index, which is in [0, n-1)
 */
inline int dsample(int n, const double *p, double u, double *F)
{
    cumsum(n, p, F);
    return get_interval(n, F, u);    
}


/**
 * Calculate the posterior
 *
 * @param n         the number of classes
 * @param lpri      the log-prior probabilities
 * @param llik      the log-likelihood with respect to each class
 * 
 * @param pos       the posterior probabilities
 */
inline void calc_posterior(int n, const double *lpri, const double *llik, 
        double *pos)
{    
    // calculate evidence (log-prior + log-lik)
    for (int i = 0; i < n; ++i)
    {
        double ev = lpri[i] + llik[i];
        pos[i] = ev;                
    }
    
    // find the maximum evidence
    double max_ev = *pos;
    for (int i = 1; i < n; ++i)
    {
        if (pos[i] > max_ev) max_ev = pos[i];
    }
    
    // shift the evidence values
    for (int i = 0; i < n; ++i)
    {
        pos[i] -= max_ev;
    }
    
    // compute exp and sum them
    double s = 0;
    for (int i = 0; i < n; ++i)
    {
        s += (pos[i] = std::exp(pos[i]));
    }
    
    // normalize to get posterior values
    double a = 1.0 / s;
    for (int i = 0; i < n; ++i)
    {
        pos[i] *= a;
    }        
}





}

#endif


